# This script shows a simple example reading and plotting the data of Fig. 7 of the paper "Measurement of the mean number of muons with energies above 500 GeV in air showers detected with the IceCube Neutrino Observatory" by the IceCube Collaboration

# imports
import numpy as np
import matplotlib.pyplot as plt

# Load the N_\mu results from data and the predictions from proton and iron simulations

E, N, Nstat, Nsysmin, Nsysplus = np.loadtxt('icecube_Nmu_sibyll21_data.txt', unpack=True)
E, N_p, N_Fe = np.loadtxt('icecube_Nmu_sibyll21_MC.txt', unpack=True)

# Make a plot similar to Fig. 7 in the paper
plt.figure()

# plot the data and uncertainties
plt.errorbar(E, N, Nstat, fmt='ok')
plt.fill_between(E, N-Nsysmin, N+Nsysplus, color='k', alpha=0.2)

# plot the MC prediction
plt.plot(E, N_p, 'r-')
plt.plot(E, N_Fe, 'b-')

# log scale
plt.xscale('log')
plt.yscale('log')

# labels
plt.xlabel('$E$ / GeV')
plt.ylabel('$N_\\mu$ ($E_\\mu > 500\,\mathrm{GeV}$)')

# indicate these are the results obtained assuming Sibyll 2.1
plt.title('Sibyll 2.1')

# save plot
figname = 'example_plot_Nmu500_sibyll21.png'
plt.savefig(figname)

print(f"Saved plot to {figname}")

